Our initial species selection includes 11 single species groups from the NOBA model. We looked at diet overlap here. Below we will double check survey census outputs against the true values, and then test some survey specifications that may be reasonable for multispecies model testing. Finally we will apply survey functions to the detailed diet data.
fgs <- load_fgs(here("atlantisoutput", "NOBA_March_2020"), "nordic_groups_v04.csv")
lname <- data.frame(Latin = c("*Hippoglossoides platessoides*",
"*Reinhardtius hippoglossoides*",
"*Scomber scombrus*",
"*Melongrammus aeglefinus*",
"*Pollachius virens*",
"*Sebastes mentella*",
"*Micromesistius poutassou*",
"*Clupea harengus*",
"*Gadus morhua*",
"*Boreogadus saida*",
"*Mallotus villosus*"),
Code = c("LRD", "GRH", "MAC", "HAD", "SAI", "RED",
"BWH", "SSH", "NCO", "PCO", "CAP")
)
sppsubset <- merge(fgs, lname, all.y = TRUE)
spptable <- sppsubset %>%
arrange(Index) %>%
select(Name, Long.Name, Latin)
knitr::kable(spptable, col.names = c("Model name", "Full name", "Latin name"))
| Model name | Full name | Latin name |
|---|---|---|
| Long_rough_dab | Long rough dab | Hippoglossoides platessoides |
| Green_halibut | Greenland halibut | Reinhardtius hippoglossoides |
| Mackerel | Mackerel | Scomber scombrus |
| Haddock | Haddock | Melongrammus aeglefinus |
| Saithe | Saithe | Pollachius virens |
| Redfish | Redfish | Sebastes mentella |
| Blue_whiting | Blue whiting | Micromesistius poutassou |
| Norwegian_ssh | Norwegian spring spawning herring | Clupea harengus |
| North_atl_cod | Northeast Atlantic cod | Gadus morhua |
| Polar_cod | Polar cod | Boreogadus saida |
| Capelin | Capelin | Mallotus villosus |
Note that survey config files must now name the survey.
NOBA2config.R looks like this (adjusted from Alfonso’s original):
d.name <- here("atlantisoutput","NOBA_march_2020")
functional.groups.file <- "nordic_groups_v04.csv"
biomass.pools.file <- "nordic_biol_v23.nc"
biol.prm.file <- "nordic_biol_incl_harv_v_007_3.prm"
box.file <- "Nordic02.bgm"
initial.conditions.file <- "nordic_biol_v23.nc"
run.prm.file <- "nordic_run_v01.xml"
scenario.name <- "nordic_runresults_01"
bioind.file <- "nordic_runresults_01BiomIndx.txt"
catch.file <- "nordic_runresults_01Catch.txt"
annage <- TRUE
fisheries.file <- "NoBAFisheries.csv"
omdimensions.R standardizes timesteps, etc. (this is part of atlantisom and should not need to be changed by the user):
#survey species inherited from omlist_ss
survspp <- omlist_ss$species_ss
# survey season and other time dimensioning parameters
# generalized timesteps all models
noutsteps <- omlist_ss$runpar$tstop/omlist_ss$runpar$outputstep
timeall <- c(0:noutsteps)
stepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutinc
midptyr <- round(median(seq(0,stepperyr)))
# model areas, subset in surveyconfig
allboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc
# survey selectivity (agecl based)
sp_age <- omlist_ss$funct.group_ss[, c("Name", "NumCohorts", "NumAgeClassSize")]
# should return all age classes fully sampled (Atlantis output is 10 age groups per spp)
n_age_classes <- sp_age$NumCohorts
# changed below for multiple species NOTE survspp alphabetical; NOT in order of fgs!!
# this gives correct names
age_classes <- sapply(n_age_classes, seq)
names(age_classes)<-sp_age$Name
n_annages <- sp_age$NumCohorts * sp_age$NumAgeClassSize
# changed below for multiple species
annages <- sapply(n_annages, seq)
names(annages)<-sp_age$Name
mssurvey_spring.R and mssurvey_fall.R configure the fishery independent surveys (in this census test, surveys sample all model polygons in all years and have efficiency of 1 for all species, with no size selectivity):
# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity
# Survey name
survey.name="BTS_spring_allbox_effic1"
#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5
#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May,
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)
survey_sample_time <- 1 # spring survey
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
total_sample, by=timestep)
survtime <- survey_sample_full
# survey area
# should return all model areas
survboxes <- allboxes
# survey efficiency (q)
# should return a perfectly efficient survey
surveffic <- data.frame(species=survspp,
efficiency=rep(1.0,length(survspp)))
# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(names(age_classes), each=n_age_classes),
# agecl=rep(c(1:n_age_classes),length(survspp)),
# selex=rep(1.0,length(survspp)*n_age_classes))
# for annage output uses names(annages) NOT alphabetical survspp
survselex <- data.frame(species=rep(names(annages), n_annages), #
agecl=unlist(sapply(n_annages,seq)),
selex=rep(1.0,sum(n_annages)))
# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))
# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1
# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200
# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity
# Survey name
survey.name="BTS_fall_allbox_effic1"
#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5
#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May,
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)
survey_sample_time <- 3 # fall survey
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
total_sample, by=timestep)
survtime <- survey_sample_full
# survey area
# should return all model areas
survboxes <- allboxes
# survey efficiency (q)
# should return a perfectly efficient survey
surveffic <- data.frame(species=survspp,
efficiency=rep(1.0,length(survspp)))
# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(survspp, each=n_age_classes),
# agecl=rep(c(1:n_age_classes),length(survspp)),
# selex=rep(1.0,length(survspp)*n_age_classes))
# for annage output
survselex <- data.frame(species=rep(names(annages), n_annages), #
agecl=unlist(sapply(n_annages,seq)),
selex=rep(1.0,sum(n_annages)))
# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))
# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1
# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200
msfishery.R configures the fishery dependent data:
# Default fishery configuration here is a census
# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss
#Number of years of data to pull
nyears <- 50
#Atlantis initialization period in years
burnin <- 30
# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc
# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample) #total_sample defined in sardinesurvey.R
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[fstepperyr], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R
fishtime <- fish_times
# fishery sampling area
# should return all model areas, this assumes you see everything that it caught
fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# effective sample size needed for sample_fish
# this effective N is divided by the number of annual timesteps below, so 200 per time
# use as input to the length samples, ages can be a subset
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))
#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr
# fishery catch cv can be used in sample_survey_biomass
# perfect observation
fish_cv <- data.frame(species=survspp, cv=rep(0.01,length(survspp)))
These survey and fishery configurations were used previously to generate saved census survey output.
Using read_savedsurvs() we get back the same objects generated by the wrapper functions om_index() and om_comps().
source(here("config", "NOBA2config.R"))
omlist_ss <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
source(here("config/omdimensions.R"))
# add this paste0(d.name, "/mscensus1") and add option to name save directory in wrappers
survObsBiom <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survB')
age_comp_data <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survAge')
len_comp_data <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survLen')
wtage <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survWtage')
annage_comp_data <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survAnnAge')
annage_wtage <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survAnnWtage')
# no function needed for saved fishery output as long as there is only one fishery specification file
catchbio_ss <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishCatch.rds")))
catchbio_ss <- catchbio_ss[[1]] #in same time unit as surveys
fish_age_comp <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]
#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- fish_age_comp$time/fstepperyr
fish_len_comp_data <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsLenComp.rds")))
#len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]
#add this to om_indices function so that this has years when read in
#fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))
fishwtage <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsWtAtAge.rds")))
fishwtage <- fishwtage[[1]] #need to change time units?
fish_annage_comp <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_annage_comp <- fish_annage_comp[[1]]
#add this to om_indices function so that this has years when read in
#fish_annage_comp$time <- fish_annage_comp$time/fstepperyr
fish_annage_wtage <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsFullWtAtAge.rds")))
fish_annage_wtage <- fish_annage_wtage[[1]] #need to change time units?
Crunch a few things to get true age comps
# age classes
totNagecl <- omlist_ss$truenums_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
totN <- totNagecl %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueNagecl <- merge(totNagecl, totN)
# annual ages
totNage <- omlist_ss$truenumsage_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
#this should be exactly the same as totN but just making sure...
totN2 <- totNage %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueNage <- merge(totNage, totN2)
# fishery catch at age class
totCagecl <- omlist_ss$truecatchnum_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
totCAAcl <- totCagecl %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueCAAcl <- merge(totCagecl, totCAAcl)
# fishery catch at annual age
totCage <- omlist_ss$truecatchage_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
totCAA <- totCage %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueCAA <- merge(totCage, totCAA)
My plotting functions: too specialized to include in the package?
# plot biomass time series facet wrapped by species
plotB <- function(dat, truedat=NULL){
ggplot() +
geom_line(data=dat, aes(x=time/stepperyr,y=atoutput, color="Survey Biomass"),
alpha = 10/10) +
{if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True B"), alpha = 3/10)} +
theme_tufte() +
theme(legend.position = "top") +
xlab("model year") +
ylab("tons") +
labs(colour=scenario.name) +
facet_wrap(~species, scales="free")
}
# make a catch series function that can be split by fleet? this doesnt
# also note different time (days) from model timestep in all other output
plotC <- function(dat, truedat=NULL){
ggplot() +
geom_line(data=dat, aes(x=time/365,y=atoutput, color="Catch biomass"),
alpha = 10/10) +
{if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True Catch"), alpha = 3/10)} +
theme_tufte() +
theme(legend.position = "top") +
xlab("model year") +
ylab("tons") +
labs(colour=scenario.name) +
facet_wrap(~species, scales="free")
}
# plot length frequencies by timestep (one species)
plotlen <- function(dat, truedat=NULL){
ggplot(dat, aes(upper.bins)) +
geom_bar(aes(weight = atoutput)) +
theme_tufte() +
xlab("length (cm)") +
ylab("number") +
labs(subtitle = paste(scenario.name,
dat$species)) +
facet_wrap(~time, ncol=6, scales="free_y")
}
# plot numbers at age by timestep (one species)
Natageplot <- function(dat, effN=1, truedat=NULL){
ggplot() +
geom_point(data=dat, aes(x=agecl, y=atoutput/effN, color="Est Comp")) +
{if(!is.null(truedat)) geom_line(data=dat, aes(x=agecl, y=numAtAge/totN, color="True Comp"))} +
theme_tufte() +
theme(legend.position = "bottom") +
xlab("age/agecl") +
ylab("number") +
labs(subtitle = paste(scenario.name,
dat$species)) +
facet_wrap(~time, ncol=6, scales="free_y")
}
# plot weight at age time series facet wrapped by species
wageplot <- function(dat, truedat=NULL){
ggplot(dat, aes(time/stepperyr, atoutput)) +
geom_line(aes(colour = factor(agecl))) +
theme_tufte() +
theme(legend.position = "bottom") +
xlab("model year") +
ylab("average individual weight (g)") +
labs(subtitle = paste0(scenario.name)) +
facet_wrap(c("species"), scales="free_y")
}
# compare N at age and C at age between standard and ANNAGE outputs
totNageplot <- function(dat, anndat){
ggplot() +
geom_line(data=dat, aes(x=time/stepperyr,y=totN, color="Tot N cohorts"),
alpha = 5/10) +
geom_line(data=anndat, aes(x=time/stepperyr,y=totN, color="Tot N annage"), alpha = 5/10) +
theme_tufte() +
theme(legend.position = "top") +
xlab("model year") +
ylab("N") +
labs(colour=scenario.name) +
facet_wrap(~species, scales="free")
}
We are using the full time series for biomass and weight at age, and a subsample of years 30-53 for compositions. In this test we see the effect of surveying in one season, with migratory species missed or only partially represented in some surveys, and with the overall surveyed biomass being at the lower or higher end of the range depending on the time of year. (Observation error is included in the survey biomass time series.) Samples for length and age assume a total of 100000 fish were randomly sampled for each species in each survey (surveffN defined in the survey config files).
# compare with true output (all timesteps)
for(s in names(survObsBiom)){
cat(" \n##### ", s," \n")
print(plotB(survObsBiom[[s]][[1]], omlist_ss$truetotbio_ss))
cat(" \n")
}
# plots survey only
# for(s in names(survObsBiom)){
# print(plotB(survObsBiom[[s]][[1]]))
# }
for(s in names(len_comp_data)){
cat(" \n##### ", s," \n")
lcompsub <- as.data.frame(len_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
group_by(species) %>%
group_map(~ plotlen(.x), keep = TRUE)
for(i in 1:length(lcompsub)) {
print(lcompsub[[i]])
}
cat(" \n")
}
for(s in names(age_comp_data)){
cat(" \n##### ", s," \n")
acompsub <- as.data.frame(age_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
group_by(species) %>%
left_join(., trueNagecl) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
cat(" \n")
}
for(s in names(annage_comp_data)){
cat(" \n##### ", s," \n")
acompsub <- as.data.frame(annage_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
group_by(species) %>%
left_join(., trueNage) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
cat(" \n")
}
for(s in names(wtage)){
cat(" \n##### ", s," \n")
print(wageplot(wtage[[s]][[1]]))
cat(" \n")
}
# annage_wtage read in above
# annage_wtage <- annage_wtage[[1]] #this still has second list component, diagnostic plot
# annage_wtage <- annage_wtage[[1]] #still in this old version but just removed from function.
for(s in names(annage_wtage)){
cat(" \n##### ", s," \n")
print(wageplot(annage_wtage[[s]][[1]][[1]]))
cat(" \n")
}
We are using the full time series for total catch and weight at age, and a subsample of years 30-35 for compositions. Fishery length and age composition is sampled every output timestep. Fishery effN for length and age sampling was 1000 fish per year (200 per time step), so this is reflected in the composition plots.
# compare with true catch, should match
plotC(catchbio_ss, omlist_ss$truecatchbio_ss)
lcompsub <- as.data.frame(fish_len_comp_data) %>% filter(time %in% c(150:175)) %>%
group_by(species) %>%
group_map(~ plotlen(.x), keep = TRUE)
for(i in 1:length(lcompsub)) {
print(lcompsub[[i]])
}
acompsub <- as.data.frame(fish_age_comp) %>% filter(time %in% c(150:175)) %>%
group_by(species) %>%
left_join(., trueCAAcl) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
acompsub <- as.data.frame(fish_annage_comp) %>% filter(time %in% c(150:175)) %>%
group_by(species) %>%
left_join(., trueCAA) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
wageplot(fishwtage)
#diagnostic plot component still in this old version but just removed from function.
wageplot(fish_annage_wtage[[1]])
mssurvey_spring_01.R and mssurvey_fall_01.R configure the fishery independent surveys. We keep spring and fall survey timing, but now we have efficiency defined for each species and a selectivity function for each species. Our surveys exclude offshore areas. These specifications are invented for demonstration purposes, not intended to realistically mimic a particular survey:
# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity
# Survey name
survey.name="BTS_spring_nearbox_qmix_selmix"
#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5
#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May,
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)
survey_sample_time <- 1 # spring survey
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
total_sample, by=timestep)
survtime <- survey_sample_full
# survey area
# should return all model areas
# survboxes <- allboxes
# shelf boxes excluding farthest east Barents
survboxes <- c(1:8, 21, 23:31, 33, 39:44, 47:49)
# survey efficiency (q)
# should return a perfectly efficient survey
# surveffic <- data.frame(species=survspp,
# efficiency=rep(1.0,length(survspp)))
# group species to apply different efficiency (all NOBA fish and shark groups)
nontrawl <- c("Sharks_other", "Pelagic_large","Mesop_fish")
pelagics <- c("Pelagic_small","Redfish_other","Mackerel","Haddock",
"Saithe","Redfish","Blue_whiting","Norwegian_ssh","Capelin")
demersals <- c("Demersals_other","Demersal_large","Flatfish_other","Skates_rays",
"Green_halibut","North_atl_cod","Polar_cod","Snow_crab")
selflats <- c("Long_rough_dab")
# define bottom trawl mixed efficiency
ef.nt <- 0.01 # for large pelagics, reef dwellers, others not in trawlable habitat
ef.pl <- 0.2 # for pelagics
ef.dm <- 0.7 # for demersals
ef.fl <- 1.1 # for selected flatfish
# bottom trawl survey efficiency specification by species group
effnontrawl <- data.frame(species=nontrawl, efficiency=rep(ef.nt,length(nontrawl)))
effpelagics <- data.frame(species=pelagics, efficiency=rep(ef.pl,length(pelagics)))
effdemersals <- data.frame(species=demersals, efficiency=rep(ef.dm,length(demersals)))
effselflats <- data.frame(species=selflats, efficiency=rep(ef.fl,length(selflats)))
efficmix <- bind_rows(effnontrawl, effpelagics, effdemersals, effselflats)
surveffic <- efficmix %>%
filter(species %in% survspp)
# survey selectivity (true age based, flat)
# for annage output uses names(annages) NOT alphabetical survspp
survselex <- data.frame(species=rep(names(annages), n_annages), #
agecl=unlist(sapply(n_annages,seq)),
selex=rep(1.0,sum(n_annages)))
# mixed selectivity: specify for annual ages 0-10 to
# then apply that curve to agecl keeping every nth selectivity index n=NumAgeClassSize
# ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]
# flat=1 for large pelagics, reef dwellers, others not in trawlable habitat
# sigmoid 0 to 1 with 0.5 inflection ~ age 3 for pelagics, reaching 1 at age 5, flat top
# sigmoid 0 to 1 with 0.5 inflection ~ age 5 for most demersals and flatfish, reaching 1 at age 7, flat top
# dome shaped 0 to 1 at agecl 6&7 for selected demersals, falling off to 0.7 by agecl 10--didn't do
sigmoid <- function(a,b,x) {
1 / (1 + exp(-a-b*x))
}
sp_age <- sp_age %>%
mutate(n_annages = NumCohorts * NumAgeClassSize)
#
# survey selectivity specification for true ages 1-10 by species group--replace for each group in surselex
selnontrawl <- data.frame(species=rep(nontrawl, each=10),
agecl=rep(c(1:10),length(nontrawl)),
selex=rep(1.0,length(nontrawl)*10))
selpelagics <- data.frame(species=rep(pelagics, each=10),
agecl=rep(c(1:10),length(pelagics)),
selex=sigmoid(5,1,seq(-10,10,length.out=10)))
seldemersals <- data.frame(species=rep(demersals, each=10),
agecl=rep(c(1:10),length(demersals)),
selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selselflats <- data.frame(species=rep(selflats, each=10),
agecl=rep(c(1:10),length(selflats)),
selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selexmix <- bind_rows(selnontrawl, selpelagics, seldemersals, selselflats)
selexmix <- selexmix %>%
filter(species %in% survspp) %>%
rename(selex10 = selex)
survselex <- merge(survselex, selexmix, all = TRUE) %>%
filter(!is.na(selex)) %>%
mutate(selex = case_when(!is.na(selex10) ~ selex10,
is.na(selex10) ~ selex)) %>%
select(-selex10)
# now apply this for agecl selectivity that matches
# and figure out how to use in wrapper!
#ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]
survselex.agecl <- survselex %>% left_join(sp_age, by=c("species"="Name")) %>%
group_by(species) %>%
filter(agecl %in% seq(unique(NumAgeClassSize),
unique(n_annages),
length.out=unique(NumCohorts))) %>%
mutate(agecl = agecl/unique(NumAgeClassSize)) %>%
select(species, agecl, selex)
# effective sample size needed for sample_fish
# this is the number of *lengths* per species that are measured on a survey
# effective N of 100000 with no age subsample matches true age comp
surveffN <- data.frame(species=survspp, effN=rep(10000, length(survspp)))
# survey index cv needed for sample_survey_xxx
# cv = 0.1
# surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
# use this constant 0 cv for testing
surv_cv_0 <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))
# define bottom trawl survey cv by group
cv.nt <- 1.0 # for large pelagics, reef dwellers, others not in trawlable habitat
cv.pl <- 0.4 # for pelagics
cv.dm <- 0.2 # for demersals
cv.fl <- 0.2 # for selected flatfish
# specify cv by species groups
surv_cv_nontrawl <- data.frame(species=nontrawl, cv=rep(cv.nt,length(nontrawl)))
surv_cv_pelagics <- data.frame(species=pelagics, cv=rep(cv.pl,length(pelagics)))
surv_cv_demersals <- data.frame(species=demersals, cv=rep(cv.dm,length(demersals)))
surv_cv_selflats <- data.frame(species=selflats, cv=rep(cv.fl,length(selflats)))
surv_cv_mix <- bind_rows(surv_cv_nontrawl, surv_cv_pelagics, surv_cv_demersals, surv_cv_selflats)
surv_cv <- surv_cv_mix %>%
filter(species %in% survspp)
# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1
# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200
# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity
# Survey name
survey.name="BTS_fall_nearbox_qmix_selmix"
#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5
#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May,
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)
survey_sample_time <- 3 # fall survey
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
total_sample, by=timestep)
survtime <- survey_sample_full
# survey area
# should return all model areas
#survboxes <- allboxes
# shelf boxes excluding farthest east Barents
survboxes <- c(1:8, 21, 23:31, 33, 39:44, 47:49)
# survey efficiency (q)
# should return a perfectly efficient survey
# surveffic <- data.frame(species=survspp,
# efficiency=rep(1.0,length(survspp)))
# group species to apply different efficiency (all NOBA fish and shark groups)
nontrawl <- c("Sharks_other", "Pelagic_large","Mesop_fish")
pelagics <- c("Pelagic_small","Redfish_other","Mackerel","Haddock",
"Saithe","Redfish","Blue_whiting","Norwegian_ssh","Capelin")
demersals <- c("Demersals_other","Demersal_large","Flatfish_other","Skates_rays",
"Green_halibut","North_atl_cod","Polar_cod","Snow_crab")
selflats <- c("Long_rough_dab")
# define bottom trawl mixed efficiency
ef.nt <- 0.01 # for large pelagics, reef dwellers, others not in trawlable habitat
ef.pl <- 0.1 # for pelagics
ef.dm <- 0.6 # for demersals
ef.fl <- 1.1 # for selected flatfish
# bottom trawl survey efficiency specification by species group
effnontrawl <- data.frame(species=nontrawl, efficiency=rep(ef.nt,length(nontrawl)))
effpelagics <- data.frame(species=pelagics, efficiency=rep(ef.pl,length(pelagics)))
effdemersals <- data.frame(species=demersals, efficiency=rep(ef.dm,length(demersals)))
effselflats <- data.frame(species=selflats, efficiency=rep(ef.fl,length(selflats)))
efficmix <- bind_rows(effnontrawl, effpelagics, effdemersals, effselflats)
surveffic <- efficmix %>%
filter(species %in% survspp)
# survey selectivity (true age based, flat)
# for annage output uses names(annages) NOT alphabetical survspp
survselex <- data.frame(species=rep(names(annages), n_annages), #
agecl=unlist(sapply(n_annages,seq)),
selex=rep(1.0,sum(n_annages)))
# mixed selectivity: specify for annual ages 0-10 to
# then apply that curve to agecl keeping every nth selectivity index n=NumAgeClassSize
# ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]
# flat=1 for large pelagics, reef dwellers, others not in trawlable habitat
# sigmoid 0 to 1 with 0.5 inflection ~ age 3 for pelagics, reaching 1 at age 5, flat top
# sigmoid 0 to 1 with 0.5 inflection ~ age 5 for most demersals and flatfish, reaching 1 at age 7, flat top
# dome shaped 0 to 1 at agecl 6&7 for selected demersals, falling off to 0.7 by agecl 10--didn't do
sigmoid <- function(a,b,x) {
1 / (1 + exp(-a-b*x))
}
sp_age <- sp_age %>%
mutate(n_annages = NumCohorts * NumAgeClassSize)
#
# survey selectivity specification for true ages 1-10 by species group--replace for each group in surselex
selnontrawl <- data.frame(species=rep(nontrawl, each=10),
agecl=rep(c(1:10),length(nontrawl)),
selex=rep(1.0,length(nontrawl)*10))
selpelagics <- data.frame(species=rep(pelagics, each=10),
agecl=rep(c(1:10),length(pelagics)),
selex=sigmoid(5,1,seq(-10,10,length.out=10)))
seldemersals <- data.frame(species=rep(demersals, each=10),
agecl=rep(c(1:10),length(demersals)),
selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selselflats <- data.frame(species=rep(selflats, each=10),
agecl=rep(c(1:10),length(selflats)),
selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selexmix <- bind_rows(selnontrawl, selpelagics, seldemersals, selselflats)
selexmix <- selexmix %>%
filter(species %in% survspp) %>%
rename(selex10 = selex)
survselex <- merge(survselex, selexmix, all = TRUE) %>%
filter(!is.na(selex)) %>%
mutate(selex = case_when(!is.na(selex10) ~ selex10,
is.na(selex10) ~ selex)) %>%
select(-selex10)
# now apply this for agecl selectivity that matches
# and figure out how to use in wrapper!
#ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]
survselex.agecl <- survselex %>% left_join(sp_age, by=c("species"="Name")) %>%
group_by(species) %>%
filter(agecl %in% seq(unique(NumAgeClassSize),
unique(n_annages),
length.out=unique(NumCohorts))) %>%
mutate(agecl = agecl/unique(NumAgeClassSize)) %>%
select(species, agecl, selex)
# effective sample size needed for sample_fish
# this is the number of *lengths* per species that are measured on a survey
# effective N of 100000 with no age subsample matches true age comp
surveffN <- data.frame(species=survspp, effN=rep(10000, length(survspp)))
# survey index cv needed for sample_survey_xxx
# cv = 0.1
# surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
# use this constant 0 cv for testing
surv_cv_0 <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))
# define bottom trawl survey cv by group
cv.nt <- 1.0 # for large pelagics, reef dwellers, others not in trawlable habitat
cv.pl <- 0.4 # for pelagics
cv.dm <- 0.2 # for demersals
cv.fl <- 0.2 # for selected flatfish
# specify cv by species groups
surv_cv_nontrawl <- data.frame(species=nontrawl, cv=rep(cv.nt,length(nontrawl)))
surv_cv_pelagics <- data.frame(species=pelagics, cv=rep(cv.pl,length(pelagics)))
surv_cv_demersals <- data.frame(species=demersals, cv=rep(cv.dm,length(demersals)))
surv_cv_selflats <- data.frame(species=selflats, cv=rep(cv.fl,length(selflats)))
surv_cv_mix <- bind_rows(surv_cv_nontrawl, surv_cv_pelagics, surv_cv_demersals, surv_cv_selflats)
surv_cv <- surv_cv_mix %>%
filter(species %in% survspp)
# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1
# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200
msfishery_01.R configures the fishery dependent data, which has a higher observation error on the total catch and does not get biological sampling from all areas:
# Default fishery configuration here is a census
fishery.name="nearbox"
# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss
#Number of years of data to pull
nyears <- 50
#Atlantis initialization period in years
burnin <- 30
# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc
# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample) #total_sample defined in survey configs
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[fstepperyr], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R
fishtime <- fish_times
# fishery sampling area (applies to comps, not total catch)
# should return all model areas, this assumes you see everything that is caught
#fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# same as survey, only getting catch comps from nearby boxes
fishboxes <- c(1:8, 21, 23:31, 33, 39:44, 47:49)
# effective sample size needed for sample_fish
# this effective N is divided by the number of annual timesteps below, so 200 per time
# use as input to the length samples, ages can be a subset
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))
#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr
# fishery catch cv can be used in sample_fishery_totcatch to
# perfect observation is cv=0
fish_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
Model boxes show where our survey works now. We are including boxes 1:8, 21, 23:31, 33, 39:44, 47:49 in these surveys, and in fishery comps sampling:
knitr::include_graphics("visualization_Rmd/images/NOBAboxes.png")
The way sampling is set up, effective N for each survey is the number of fish measured for length, age, and average weight at age. We could introduce further realism by using the sample_ages() function to take a subsample for age composition and optionally apply an ageing error matrix. This would require changing the om_comps() wrapper: see atlantisom issue.
#NOBAom <- om_init(here("config/NOBA2config.R"))
#NOBAom_ms <- om_species(sppsubset$Name, NOBAom)
NOBAom_ms <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
# new wrapper for multiple surveys saves survey output separately
# but returns a list with all surveys, not ideal; make consistent
NOBAom_ms_ind <- om_index(usersurvey = c(here("config/mssurvey_spring_01.R"),
here("config/mssurvey_fall_01.R")),
userfishery = here("config/msfishery_01.R"),
omlist_ss = NOBAom_ms,
n_reps = 1,
save = TRUE)
NOBAom_ms_comp <- om_comps(usersurvey = c(here("config/mssurvey_spring_01.R"),
here("config/mssurvey_fall_01.R")),
userfishery = here("config/msfishery_01.R"),
omlist_ss = NOBAom_ms,
n_reps = 1,
save = TRUE)
NOTE! Need to have selectivities defined both for annual ages and for age classes. Worked in census because flat selectivity was defined for full set of annual ages and it just didn’t use past age 10 for age class output. They were all 1 anyway so it didnt matter. The updated survey config file first specifies selectivity by true age and then backs out an approximate equivalent for age classes. We will use true age classes for the real datasets but wanted to check both outputs here. om_comps() has now changed to take survselex.agecl as input to the age class comps and survselex as input to the full comps.
Both om_index() and om_comps() wrappers now name fishery output based on a fishery.name in the fishery config file.